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ABSTRACT 

Recent experimental results on the formation of molecular hydrogen on astrophysi- 
cally relevant surfaces under conditions similar to those encountered in the interstellar 
medium provided useful quantitative information about these processes. Rate equation 
analysis of experiments on olivine and amorphous carbon surfaces provided the acti- 
vation energy barriers for the diffusion and desorption processes relevant to hydrogen 
recombination on these surfaces. However, the suitability of rate equations for the sim- 
ulation of hydrogen recombination on interstellar grains, where there might be very few 
atoms on a grain at any given time, has been questioned. To resolve this problem, we 
introduce a master equation that takes into account both the discrete nature of the H 
atoms and the fluctuations in the number of atoms on a grain. The hydrogen recombi- 
nation rate on microscopic grains, as a function of grain size and temperature, is then 
calculated using the master equation. The results are compared to those obtained from 
the rate equations and the conditions under which the master equation is required are 
identified. 

Subject headings: dust — ISM; abundances — ISM; molecules — molecular processes 



1. Introduction 

The formation of molecular hydrogen in the interstellar medium (ISM) is a process of funda- 
mental importance (Duley & Williams 1984; Williams 1998). It was recognized long ago (Gould 
& Salpeter 1963) that H2 cannot form in the gas phase efficiently enough to account for its abun- 
dance. It was thus proposed that dust grains act as catalysts, where an H atom approaching the 
surface of a grain has a probability £ to become adsorbed. The adsorbed H atom (adatom) spends 
an average time £h (residence time) before leaving the surface. If during the residence time the 
H adatom encounters another H adatom, an H2 molecule will form with a certain probability. 
Various aspects of this process were addressed in extensive theoretical studies (Gould & Salpeter 
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1963; Williams 1968; Hollenbach & Salpeter 1970, 1971; Hollenbach et al. 1971; Smoluchowski 1981, 
1983; Aronowitz & Chang 1985; Duley & Williams 1986; Pirronello & Averna 1988; Sandford & 
Allamandolla 1993; Takahashi et al. 1999; Farebrother et al. 1999). In particular, Hollenbach et al. 
calculated the sticking and mobility of H atoms on grain surfaces. They concluded that tunneling 
between adsorption sites, even at temperature as low as T = 10K, provides the required mobility. 
The steady state production rate of molecular hydrogen per unit volume was expressed according 
to (Hollenbach et al. 1971) 

^H 2 = -puvuvipg, (1) 

where pn and t>H are the number density and the speed of H atoms in the gas phase, respectively, 
a is the average cross-sectional area of a grain and p g is the number density of dust grains. The 
parameter 7 is the fraction of H atoms striking the grain that eventually form a molecule, namely 
7 = £77, where r\ is the probability that an H adatom on the surface will recombine with another H 
atom to form H2. 

Recently, a series of experiments were conducted to measure hydrogen recombination in an 
ultra high-vacuum (UHV) chamber by irradiating the sample with two beams of H and D atoms 
and monitoring the HD production rate (Pirronello et al. 1997a, b, 1999). Two different substrates 
have been used: a natural olivine (a polycrystalline silicate containing Mg2Si04 and Fe2SiC"4) slab 
and an amorphous carbon sample. The substrate temperatures during hydrogen irradiation were 
in the range between 5 K and 15 K. The HD formation rate was measured using a quadrupole mass 
spectrometer both during irradiation and in a subsequent temperature programmed desorption 
(TPD) experiment in which the sample temperature was quickly ramped to over 30 K to desorb 
all weakly adsorbed species. It was found that H and D atoms adsorbed on the surface at the 
lowest irradiation temperature of 5 K form molecules during TPD only above 9 K in the case 
of olivine and above 14 K in the case of amorphous carbon. This indicates that tunneling alone 
does not provide enough mobility to H adatoms to enable recombination, and thermal activation is 
required. The experimental results were analyzed using a rate equation model (Katz et al. 1999). 
In this analysis the parameters of the rate equations were fitted to the experimental TPD curves. 
These parameters are the activation energy barriers for atomic hydrogen diffusion and desorption, 
the barrier for molecular hydrogen desorption and the probability of spontaneous desorption of 
a hydrogen molecule upon recombination. Using the values of the parameters that fit best the 
experimental results, the efficiency of hydrogen recombination on the olivine and amorphous carbon 
surfaces was calculated for interstellar conditions using the same rate equation model. It was 
found that for both samples the recombination efficiency is strongly dependent on temperature and 
exhibits a narrow window of high recombination efficiency along the temperature axis. 

It was recently pointed out that since hydrogen recombination in the interstellar space takes 
place on small grains, rate equations have a limited range of validity (Tielens 1995; Charnley et al. 
1997; Caselli et al. 1998; Shalabiea et al. 1998). This is due to the fact that these equations take 
into account only average concentrations and ignore fluctuations as well as the discrete nature of 
the H atoms. These properties become significant in the limit of very small grains and low incoming 
flux of H atoms, exactly the conditions encountered in diffuse interstellar clouds where hydrogen 
recombination on silicate and carbon surfaces is expected to be relevant. As the number of H atoms 
on a grain fluctuates in the range of 0, 1 or 2, the recombination rate cannot be obtained from the 
average number alone. This can be easily understood, since the recombination process requires at 
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least two H atoms simultaneously on the surface. Comparisons with Monte Carlo simulations have 
shown that the rate equations tend to over-estimate the recombination rate. A modified set of rate 
equations which exhibits better agreement with Monte Carlo simulations was introduced by Caselli 
et al. (1998) and applied by Shalabiea et al. (1998) to a variety of chemical reactions. In these 
equations the rate coefficients are modified in a semi-empirical way to take into account the effect 
of the finite grain size on the recombination process. 

In this paper we introduce a master equation that is particularly suitable for the simulation of 
chemical reactions on microscopic grains. It takes into account both the discrete nature of the H 
atoms as well as the fluctuations. Its dynamical variables are the probabilities Ph(Nu) that there 
are Nn atoms on the grain at time t. The time derivatives Pn(N}i), jVh = 0, 1, 2, . . . are expressed 
in terms of the adsorption, reaction and desorption terms. The master equation provides the time 
evolution of Ph(Nyi), = 0,1,2, .. ., from which the recombination rate can be calculated. We 
use it in conjunction with the surface parameters obtained experimentally, to explore the hydrogen 
recombination process on microscopic grains for grain sizes, flux and surface temperatures pertinent 
to the conditions in the interstellar medium. 

The paper is organized as follows. The rate equation model is described in Sec. 2. The master 
equation is introduced in Sec. 3. Computer simulations and results for hydrogen recombination on 
microscopic grains under interstellar conditions are presented in Sec. 4. The case of more complex 
reactions involving multiple species is considered in Sec. 5 and a summary in Sec. 6. 



2. Rate Equations for H2 Formation on Macroscopic Surfaces 

Consider an experiment in which a flux of H atoms is irradiated on the surface. If the tem- 
perature is not too low H atoms that stick to the surface perform hops as random walkers and 
recombine when they encounter one another. Let nn(t) (in monolayers [ML]) be the coverage of H 
atoms on the surface and nu 2 (t) (also in ML) the coverage of H2 molecules at time t. We obtain 
the following set of rate equations: 

= /h • (1 - «h - «h 2 ) - Wnnn - 2a H n H (2a) 

— ^— = /j,ann H -Wu 2 nu 2 . (2b) 

The first term on the right hand side of Eq. (2a) represents the flux of H atoms multiplied by the 
Langmuir-Hinshelwood rejection term. In this scheme H atoms deposited on top of H atoms or H2 
molecules already on the surface are rejected. The parameter /h represents the effective flux of 
atoms (in units of MLs~ l ), namely, the (temperature dependent) sticking coefficient £(T) of the 
bare surface is absorbed into /h- The second term in Eq. (2a) represents the desorption of H atoms 
from the surface. The desorption coefficient is 

W }1 = v-exp{-E l /k B T) (3) 

where v is the attempt rate (standardly taken to be 10 12 s _1 ), E\ is the activation energy barrier 
for desorption of an H atom and T is the temperature. The third term in Eq. (2a) accounts for the 
depletion of the H population on the surface due to recombination into H2 molecules, where 



a H = v ■ eyLp(-E /k B T) 



(4) 
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is the hopping rate of H atoms on the surface and Eq is the activation energy barrier for H dif- 
fusion. Here we assume that diffusion occurs only by thermal hopping, in agreement with recent 
experimental results (Katz et al. 1999). We also assume that there is no barrier for recombination. 
The first term on the right hand side of Eq. (2b) represents the creation of H2 molecules. The 
factor 2 in the third term of Eq. (2a) does not appear here since it takes two H atoms to form one 
molecule. The parameter \i represents the fraction of H2 molecules that remains on the surface 
upon formation, while a fraction of (1 — /z) is spontaneously desorbed due to the excess energy 
released in the recombination process. The second term in Eq. (2b) describes the desorption of H2 
molecules. The desorption coefficient is 

W H2 =v-exp(-E 2 /k B T), (5) 

where E2 is the activation energy barrier for H2 desorption. The H2 production rate rn 2 (ML s _1 ) 
is given by: 

r H2 = (1 - n) ■ a H nn + Wft 2 rift 2 . (6) 

This model can be considered as a generalization of the Polanyi-Wigner equation [see e.g. Chan et 
al. (1978)]. It provides a description of both first order and second order desorption kinetics for 
different regimes of temperature and flux. 

The model described by Eqs. (2) was used by Katz et al. (1999) to analyze the results of the 
TPD experiments (Pirronello et al. 1997a,b, 1999). The values of the parameters Eq, E\, E2, and 
/j,, that best fit the experimental results were obtained. For the olivine sample it was found that 
Eq = 24.7 meV, E\ = 32.1 meV, Ei = 27.1 meV and fi = 0.33, while for the amorphous carbon 
sample Eq = 44.0 meV, E x = 56.7 meV, E 2 = 46.7 meV and /x = 0.413. 

The model [Eqs. (2)] was then used in order to describe the steady state conditions that 
are reached when both the flux and the temperature are fixed. The steady state solution is then 
easily obtained by setting drift /dt = and dnft 2 /dt = and solving the quadratic equation for nft 
(Biham et al. 1998; Katz et al. 1999). In case that the Langmuir-Hinshelwood rejection term can 
be neglected, the steady-state coverages are 



-Wft + JW£:+8aftfft 
U « = -4a~ft (7a) 

nft 2 = (wi + Aciftfft -Wft^W u + Saftfft^j . (7b) 

More complicated expressions are obtained when the rejection term is taken into account (Katz et 
al. 1999). The recombination efficiency 77 is defined as the fraction of the adsorbed H atoms that 
desorb in the form of H2 molecules, namely 

Note that under steady state conditions 77 is limited to the range < 77 < 1. 

By varying the temperature and flux over the astrophysically relevant range the domain in 
which there is non-negligible recombination efficiency was identified. It was found that the recom- 
bination efficiency is highly temperature dependent. For each of the two samples there is a narrow 
window of high efficiency along the temperature axis, which shifts to higher temperatures as the 
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flux is increased. For the astrophysically relevant flux range the high efficiency temperature range 
for olivine was found to be between 7 — 9K, while for amorphous carbon it is between 12 — 16 K. 

Note that in steady state the dependence of the production rate rn 2 on \i and Wn 2 is only 
through the Langmuir-Hinshelwood rejection term. This is easy to understand since the parameter 
\x only determines what fraction of the H2 will desorb upon formation and what fraction will desorb 
thermally at a later time. The desorption rate Wh 2 determines the coverage of H2 molecules at 
steady state. As long as the coverage of H and H2 is low, the Langmuir-Hinshelwood rejection term 
is small and fi and Wn 2 have little effect on the production rate rn 2 • Under interstellar conditions 
the coverage is expected to be low. Therefore, the master equation presented below, that is required 
only when the number of atoms on the grain is small, does not include the rejection term. Note, 
however, that at lower temperatures in which H atoms are immobile (thus recombination and 
desorption are suppressed) they may accumulate on the surface and reach a high coverage. 



3. Master Equation for H2 Formation on Small Grains 

We will now consider the formation of H2 molecules on small dust grains in interstellar clouds. 
In this case it is more convenient to rescale our parameters such that instead of using quantities per 
unit area - the total amount per grain will be used. The number of H atoms on the grain is denoted 
by Ah- Its expectation value is given by (Ah) = S • «h where S is the number of adsorption sites 
on the grain. Similarly, the number of H2 molecules on the grain is Ah 2 and its expectation value 
is (Ah 2 ) = S ■ ri}i 2 (we assume that each adsorption site can adsorb either an H atoms or an H2 
molecule). The incoming flux of H atoms onto the grain surface is given by Fh = S ■ fu (atoms 
s _1 ). The desorption rates Wr and Wu 2 remain unchanged. The hopping rate an (hops s _1 ) is 
replaced by An = au/S which is approximately the inverse of the time t s required for an atom 
to visit nearly all the adsorption sites on the grain surface. This is due to the fact that in two 
dimensions the number of distinct sites visited by a random walker is linearly proportional to the 
number of steps, up to a logarithmic correction (Montroll &; Weiss 1965). The H2 production rate 
of the single grain is given by Ru 2 = S • rn 2 (molecules s -1 ). The rate equations (neglecting the 
Langmuir-Hinshelwood rejection term) will thus take the form 

' /(All> = /•„ Hh .Yh. 2.1„ A,, 2 (9a) 



dt 
d(N H2 ) 
dt 



Fh 2 + Mh(A h ) 2 - WuANn,), (9b) 



where the first term in (9b) accounts for the flux of hydrogen molecules from the gas phase that 
are adsorbed on the grain surface. While for large grains Eqs. (9) provide a good description of 
the recombination processes, in the limit in which the number of atoms on the grain becomes small 
they are not suitable anymore. 

In order to to resolve this problem we will now introduce a different approach based on a 
master equation that is suitable for the study of H2 formation on small grains. Each grain is 
exposed to a flux Fh of H atoms. At any given time the number of H atoms adsorbed on the grain 
may be Ah = 0, 1, 2, — The probability that there are Ah hydrogen atoms on the grain is given 
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by Ph{N h ), where 



oc 



2 Ph(Nk) = 1. 



(10) 



AT H =0 



The time derivatives of these probabilities, Ph(A/h), are given by 



Ai(0) = -P H Ph(0) + W h Pr(1) + 2 ■ 1 • A K P K (2) 

P H (1) = F h [Ph(0)-Ph(1)]+^h[2Ph(2)-Ph(1)]+3-2-A h Ph(3) 

Ph(2) = F h [Ph(1)-Ph(2)]+^h[3Ph(3)-2P h (2)] 

+ A h [4-3-Ph(4)-2-1-Ph(2)] 

P H (N H ) = F H [p H (N H -l)-P H (N H )]+W il [(N H + l)P ll (N u + l)-N ll P ll (N u )} 

+ A H [(N H + 2){N K + l)P H {N H + 2)-N H (N K -l)P H {N H )}. (11) 



Each of these equations includes three terms. The first term describes the effect of the incoming 
flux Fh on the probabilities. The probability Ph(A/h) increases when an H atom is adsorbed on a 
grain that already has Ah — 1 adsorbed atoms [at a rate of PhPh(A"h — 1)], and decreases when 
a new atom is adsorbed on a grain with iVn atoms on it [at a rate of PhPh(-^h)]- The second 
term includes the effect of desorption. An H atom desorbed from a grain with Nn adsorbed atoms 
decreases the probability Ph(Ah) [at a rate of NhWhPh(Nh), where the factor A/h is due to the 
fact that each of the A/h atoms can desorb] and increases the probability Ph(A ? h — 1) at the same 
rate. The third term describes the effect of recombination on the number of adsrobed H atoms. 
The production of one molecule reduces this number from A/h to A/h — 2. For one pair of H atoms 
the recombination rate is proportional to the sweeping rate An multiplied by 2 since both atoms 
are mobile simultaneously. This rate is multiplied by the number of possible pairs of atoms, namely 
A/h (A/h — l)/2. Note that the equations for Ph(0) and Ph(1) do not include all the terms, because 
at least one H atom is required for desorption to occur and at least two for recombination. The 
rate of formation of H2 molecules, Ph 2 (molecules s _1 ), is thus given by 



N H =2 

Note that the sum could start from A^h = since the first two terms vanish. The recombination 
efficiency is given by 



The probability that there are Ah 2 hydrogen molecules on the grain is given by Ph 2 (A^h 2 ). The 
time evolution of these probabilities is given by 



P H2 (0) = -Ph 2 Ph 2 (0) + ^h 2 Ph 2 (1)-mP H2 Ph 2 (0) 

Ph 2 (1) = Ph 2 [Ph 2 (0) - Ph 2 (1)] + Wu 2 [2P H2 (2) - Ph 2 (1)] + /uPh 2 [Ph 2 (0) - Ph 2 (1)] 
P H2 (2) = P H2 [Ph 2 (1) - Ph 2 (2)] + W n2 [3P Ha (3) - 2P Ha (2)] + ^ [Ph 2 (1) - Ph 2 (2)] 



00 




(12) 
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Ph 2 (Nh 2 ) = F H2 [P H2 m 2 - 1) - Pu 2 (Nu 2 )} + W H , [(N H2 + l)P H2 (iVH 2 + 1) - N H2 Pa 2 (N H2 )] 
+ nRH a [Pa a i.N Ha -l)-Pk a (NK 2 )] (14) 

where -Fh 2 (molecules s _1 ) is the flux of H2 molecules that stick on the grain surface, Wh 2 is the 
desorption rate of molecules from the surface (which is inversely proportional to their residence 
time, namely in 2 = 1/Wh 2 )- Each of these equations includes three terms, describing the effects of 
an incoming H2 flux, desorption and recombination, respectively. 

Note that the fact that we ignored the Langmuir-Hinshelwood rejection term allowed us to split 
the master equation into two parts: Eq. (11) for the H atoms and Eq. (14) for the H2 molecules. 
Moreover, Eq. (11) does not depend on the distribution of Nn 2 , while Eq. (14) depends only on the 
first and second moments of the distribution of Ah- The most general case, in which the rejection 
term is included, would require to use a master equation for the joint probability distribution 
i r H&H 2 (-^H ) Ah 2 ), which is clearly much more complicated. 

The expectation value for the number of H atoms on the grain is 

00 

(N R )= Yl ^hPh^h) (15) 

N H =0 

and the expectation value for the number of molecules is 

00 

<^h 2 ) = N H2 Ph 2 (N H2 ). (16) 

^Vh 2 =0 

The time dependence of these expectation values, obtained from Eqs. (11) and (14), is given by 

d(N H ) 



dt 

d{N U2 ) 



= Fn-WniNn) - 2A K {N U (N U - 1)) (17) 
= F H2 +fiA H (N H (N H -l))-W H2 {N H2 ), (18) 



dt 

and the recombination rate i?H 2 (molecules s^ 1 ) for the grain is 

Rn 2 = (1 - »)A H (N H (N H - 1)) + W Ha (N Ha ). (19) 

These equations resemble the rate equations (9) apart from one important difference: the recombi- 
nation term (-/Vh) 2 is replaced by (-/Vh 2 ) — (-Vh)- On a macroscopically large grain it is expected that 
the difference between these two terms will be small and Eqs. (9) would provide a good description 
of the system. However, on a small grain, where (-/Vh) is small these two terms are significantly 
different and it is necessary to use the master equation rather than the rate equations. 

In principle the master equation consists of infinitely many equations for each atomic or molec- 
ular specie. In practice, for each specie such as atomic hydrogen we simulate a finite number of 
equations for Fh(Ah), Ah = 1, ... , AW, where Pr(Nb_) = for Nu > AW- Obviously, AW can- 
not exceed the number of adsorption sites, S, on the grain. In the equations for Ph (N max — 1) and 
for Pf/(A r max ), the terms that couple them to Pjj (JV max + 1) and PH(N m & x + 2) are removed (these 
terms describe the flow of probability from Pn(Nft), Ah > N m£LX to -Ph(Ah)> Ah < A^ max ). Terms 
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such as FhPii(Nii) that describe probability flow in the opposite direction are also removed. The 
latter terms are evaluated separately and frequently during the integration of the master equation 
in order to examine whether iV max should be increased. The condition for adding more equations 
is typically Pn(N max + 1) > e at a certain time t, where e is a small parameter, suitably chosen 
according to the desired precision. 

Note that the master equation is typically needed when (Ah) is of order unity. Under such 
conditions most of the probability Ph(Ah) is concentrated at small values of Ah and therefore 
Amax is expected to be small. In a time dependent simulation, when (Ah) increases reaching the 
limit (Nu) 3> 1 (thus requiring iV max 3> 1) the master equation can be easily replaced by the rate 
equations, during the run. One simply has to evaluate (Nu) at a certain time t and from that 
point to continue the run using the rate equations. The opposite move of switching from the rate 
equations to the master equation (when (Nu) decreases towards (Nu) « 1) is nearly as simple. 
One has to pick as an initial condition for the master equation a narrow distribution Pn(Nft) 
that satisfies the average (Ah) given by the rate equations, and after some relaxation time it will 
converge to the proper distribution. In simulations of more complex reactions involving multiple 
species, the coupling between different species typically involves only averages such as (Ah) (this 
is an approximation that will be discussed below). Therefore, one can simultaneously use the rate 
equations for some species and the master equation for others, according to the criteria mentioned 
above. 

To couple the master equation, consisting of Eqs. (11) and (14), to the densities in the gas 
phase we will consider the densities pu (atoms cm' 3 ) of H atoms and pu 2 (molecules cm -3 ) of H2 
molecules in the gas phase. The incoming fluxes onto the surface of a single grain can be expressed 
as Fh = Phvhc and Fh 2 = Ph 2 v h 2 °' where vu 2 is the average speed of an H2 molecule in the gas 
phase. The time derivatives of the densities are given by 

pn = [-F n + W n (N u )]p g (20) 
Pu 2 = [-F H2 + (l-p)A H (N H (N H -l)) + W H2 (N H2 )}p g . (21) 

4. Simulations and Results 

To examine the effect of the finite grain size on the recombination rate of hydrogen in the 
interstellar medium we performed simulations of the recombination process using the master equa- 
tion [that consists of Eqs. (11) and (14)] and compared the results to those obtained from the rate 
equations (9). Since we focus here on steady state conditions, only the part of the master equation 
included in Eq. (11) needs to be integrated, and the recombination rate is given by Eq. (12). For 
non-steady state conditions, Eq. (14) would be essential in order to evaluate the time dependent 
recombination rate. The parameters we have used are given below. Assuming, for simplicity, a 
spherical grain of diameter d we obtain a cross section of a = ird 2 /4. The estimate of the number 
of adsorption sites on the grain was based on the experimental data for the olivine and amorphous 
carbon surfaces (Pirronello et al. 1997a, 1997b, 1999) using the following procedure. The flux of 
the H and D beams was estimated as b = 10 12 (atoms cm" 2 s" 1 ). The beams passed through a 
chopper that reduced their flux by a factor of c = 20. A measurement of the flux in units of ML 
per second was done using the data for the total HD yield vs. exposure time [Fig. 3 in Pirronello 
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et al. (1997a)]. The theoretical Langmuir-Hinshelwood mechanism provides a prediction for the 
coverage of adsorbed atoms after irradiation time t, which is 

n H (t) = 1 - exp(-/ H ■ t). (22) 

Fitting the total HD yield to this expression we obtained good fits that provide the flux values 
fu = 2.7 • 1CP 4 (in ML s _1 ) for the olivine experiment and fu = 9.87 • 1CP 4 for the amorphous 
carbon experiment. From these two measurements we obtain the density of adsorption sites (sites 
cm~ 2 ) 

s = * (23) 

For the olivine surface it is found that s = 2 ■ 10 14 and for the amorphous carbon surface s = 5 ■ 10 13 
(sites cm" 2 ). 

Observations indicate that there is a broad distribution of grain sizes, that roughly resembles 
power-law behavior, in the range of 10 -6 cm < d < 10~ 4 cm (Mathis et al. 1977, 1996; O'Donnell 
& Mathis 1997). The number of adsorption sites on a (spherical) grain is given by S = ird 2 s. In 
the simulations we focus on diffuse clouds and use as a typical value for the density of H atoms 
pu = 10 (atoms cm" 3 ). The temperature of the H gas is taken as T = 100K. The typical velocity 
of H atoms in the gas phase is given by (see e.g. Landau & Lifshitz 1980) 



/ 8 kBT (OA\ 

VH = \/ (24) 

V 7r m 

where m = 1.67 • 10" 24 ( gram) is the mass of an H atom. We thus obtain vh = 1-45 • 10 5 (cm s 1 ). 
The density of grains is typically taken as p g = 10~ 12 /9h and hence in our case p g = 10~ n (grains 
cm" 3 ). The sticking probability of H atoms onto the grain surface is taken as £ = 1. Experimental 
results indicate that the sticking probability is close to 1 for temperatures below about 10K and 
possibly somewhat lower at higher temperatures. Since there is no high quality experimental data 
for the temperature dependence £(T), and in order to simplify the analysis we chose £ = 1. 

We will now analyze the processes that take place on a single grain using numerical integration 
of the master equation and comparison to the rate equations. The flux of H atoms onto the grain 
surface is given by fu = puvu/{^s) (ML s _1 ), where the factor of 4 in the denominator is the ratio 
between the surface area and the cross section for a spherical grain. Using the parameters above 
we obtain that fu = 0.18 • 10" 8 for olivine and fu = 0.73 • 10" 8 for amorphous carbon. The total 
flux on a grain of diameter d is given by Fh = fu ■ S (atoms s" 1 ). 

In Fig 1 we present the recombination efficiency rj for an olivine grain, under steady state 
conditions, as a function of the grain temperature. The solid line shows the results obtained 
from the rate equations, showing a range of very high efficiency between 7 — 9 K and a tail of 
decreasing efficiency between 9 — 10 K. The results of the master equation are shown for grains of 
diameters d = 10" 5 cm (O) and d = 10~ 6 cm (x). In this case the total flux on a grain amounts 
to Fh = 10~ 3 and 1CT 5 (atoms s _1 ) for the larger and smaller grain sizes, respectively. It is 
found that for the larger grain there is good agreement between the master equation and rate 
equation results. However, for the smaller grain the master equation predicts significantly lower 
recombination efficiency for temperatures higher than 8 K. Note that in order to produce the rate 
equation results for the entire temperature range of 5 - 15 K we had to include in these equations 
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the Langmuir-Hinshelwood rejection term (Katz et al. 1999). In the low temperature limit, where 
atoms are immobile (and the recombination rate decreases to zero), this is necessary in order to 
make sure that the coverage does not exceed 1 ML. For higher temperatures, where the comparison 
with the master equation results is made, the coverage is low and the rejection term makes no 
difference. From the results for the recombination efficiency one obtains the production rate of H 2 
molecules 7Zh 2 = ^FuPgiJ (molecules cm' 3 s _1 ) released to the gas phase by a density p g of grains. 

In Fig. 2 we present the expectation value for the number (Nh) of H atoms on an olivine 
grain as a function of the grain diameter d. The results (O) ar e obtained from the master equation 
under steady state conditions at T = 10K. The solid line is simply a guide to the eye. As may be 
expected, (Nh) ~ d 2 , namely, it is proportional to the surface area of the grain. It is observed that 
for a grain diameter smaller than about 10 -5 (cm) the expectation value (Nh) decreases below 
one H atom on the grain. Under these conditions significant deviations are expected between the 
recombination efficiency predicted by the rate equations and the correct results obtained from the 
master equation. 

The recombination efficiency rj for an olivine grain as a function of the grain size (•) is shown 
in Fig. 3. The temperature and the flux are identical to those used in Fig. 2. The dashed line shows 
the recombination efficiency obtained from the rate equations under similar conditions, which is 
independent of the grain size. It is observed that for grain diameter smaller than about 10 -5 
(cm) the recombination efficiency sharply drops below the rate-equation value. This is due to the 
fact that in this range (Nh) < 2, hence most often an H atom resides alone on the grain and no 
recombination is possible. 

In Fig. 4 we present the distribution Pu(Nn) on olivine grains of diameters d = 10~ 5 cm (O) 
and d = 1CT 6 cm (x). The results were obtained from the master equation under steady state 
conditions at T = 9 K. For the larger grain the distribution exhibits a peak around (Nh) — 14 and 
the rate equations are expected to apply. However, for the smaller grain the highest probability is 
for having no H atoms at all on the grain and (Nh) < 1- Under these conditions the rate equations 
are expected to highly over-estimate the recombination efficiency. Indeed, this can be observed in 
Fig. 1, in which the recombination efficiency for the smaller grain size at 9 K is much lower than 
the rate equation result. 

The results for the recombination efficiency on amorphous carbon grains are qualitatively 
similar to those shown here for olivine. The temperature range of very high efficiency is between 
12 — 16 K and the tail of decreasing efficiency is between 16 — 18 K [see Figs. 6 and 7 in Katz 
et al. (1999)]. This high efficiency window is within the relevant temperature range for diffuse 
cloud environments. It is observed that for temperatures higher than about 15 K and grain sizes 
smaller than 10~ 5 (cm) the master equation predicts lower recombination efficiency than the rate 
equations. 

5. More Complex Reactions of Multiple Species 

The master equation introduced above can be extended to describe more complex situations 
that involve chemical reactions with multiple species. We chose to demonstrate this procedure for 
the reactions involving oxygen and hydrogen on dust grains, previously studied by Caselli et al. 



- 11 - 



P T T Wtt//Vtt\ 9.4tt/7Vtt\ 2 








F - (A H + A )(N U )(N ) - 2A (N ) 2 


(25b) 


fiA H (N H ) 2 -W H2 (N H2 ). 


(25c) 


(A R + A )(N U }(N ) - A H {N H ){N OH ) 


(25d) 


Ao(N ) 2 


(25e) 


Au(N n)(N H } 


(25f) 



(1998). The rate equations that describe these reactions are 

d(N u ) 
dt 

d(Np) 

dt 
d{N U2 ) 

dt 
d(N ou ) 

dt 

d(Np 2 } 
dt 

d(N H2 o) 
dt 

where No is the number of oxygen atoms on the grain and Ao is their sweeping rate. The numbers 
of H2, OH, O2 and H2O molecules on the grain are given by Nn 2 , Nou, N<j 2 and Nr 2 o, respectively. 
The flux of O atoms adsorbed on the grain surface is given by Fo (atoms s _1 ) while, for simplicity, 
the adsorption of the four molecular species from the gas phase is neglected. Since the oxygen 
atoms and the molecules they form are chemisorbed on the surface they are unlikely to desorb 
for the surface temperatures considered here. Their desorption coefficients are thus neglected. 
It is also assumed that the diffusion of the four molecular species is negligible. Except for the 
OH molecule, this assumption is inconsequential since the three other molecular species do not 
participate in any subsequent reactions. Furthermore, it is typically assumed that the diffusion 
rate of oxygen atoms is much slower than of hydrogen, namely Ao <C ^h- In the analysis below 
we will assume that Ao = 0, namely the reaction between H and O is driven only by the diffusion 
rate of the hydrogen atoms. Under this assumption and for low coverage of O atoms on the grain, 
the production of oxygen molecules is suppressed and Eq. (25e) can be ignored. The assumption 
that H atoms are much more mobile than heavier atomic species such as C, O and N implies that 
the diffusion of H atoms is dominant in other reactions that take place on grain surfaces besides 
hydrogen recombination. Thus, the activation energies obtained by Katz et al. (1999) may be used 
not only for the hydrogen recombination process but for a large number of other reactions on dust 
grain surfaces. However, one should carefully examine the possibility that some of these reactions 
involve an additional activation energy associated with the reaction itself. 

For simplicity we will now neglect the reaction between H atoms and OH molecules which 
generates H2O molecules. We will also assume, for simplicity, that H2 molecules desorb into the 
gas phase immediately upon formation, namely that fi = 0. Under these conditions the rate 
equations will simplify into 



d(N R ) 

dt 
d(Np) 

dt 
d(N OH ) 
dt 



= F R -W R (N U )-2A R (N R ) 2 -A U (N R }(N ) (26a) 
= F -A R (N H )(N ) (26b) 
= A K (N U )(N ). (26c) 



The part of the master equation describing the time evolution of the population of H atoms on the 
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jrain surface takes the form 

P H (0) = -FhPh(O) + W H P H (1) + 2 • 1 • A H P H (2) + A K Pn(l){N ) 

P H (1) = F h [Ph(0)-P h (1)] + Wh[2Ph(2)-Ph(1)] 

+ 3 • 2 • A H P H (3) + A H [2P H (2) - P H (1)] (iV ) 

Ph(2) = P h [Ph(1) -Ph(2)] +VFh[3Ph(3) -2Ph(2)] 

+ A H [4 ■ 3 ■ P H (4) - 2 • 1 • Ph(2)] + An [3P H (3) - 2P H (2)] (N ) 

Pn(Nu) = F U [PuiNn - 1) - Ph(7V H )] + VFh [(JV h + l)fb(JV H + 1) - A^hPh^h)] 

+ A H [(JV H + 2)(iV H + l)P H (A f H + 2) - A^h - l)P H (A r H )] 

+ A H [(JVh + l)fk(^ H + 1) - N u Pa(N u )] (N Q ) (27) 



where 

oo 

(A^o) = Yl N oPo(No). (28) 

N =l 

The equations for oxygen atoms are 

P o (0) = -F o P H (0) + A R P (1)(N R ) 

Po(l) = F [P o (0) - Po(l)] + A u [2P g (2) - P G (1)] (JV H > 

Po(2) = Po[Po(l)-Po(2)]+A H [3Po(3)-2Po(2)](iVH) (29) 

Po(N ) = F [Po(No - 1) - Poi.No)] + A n [{N Q + 1)P (N + 1) - N Po(No)} (N u ). 



where 

oo 

(Nh) = Yl N R P H (N U ). (30) 

JV H =1 

The equations describing the distribution of the number of OH molecules are 
Poh(0) = -A R P OR (0)(N o )(N R ) 

P OH (l) = +A H [PoH(0)-P O H(l)](iVo}(iVH) 

P OH (2) = +A k [P k(1)-Pok(2)](No)(N k ) (31) 
Pou(Nou) = +A n [P r(N r - 1) - Pqh^oh)] (Nq) (Nu). 



The rate of formation of H2 molecules is given by Eq. (12). The rate of formation of OH molecules 
is given by Poh = Ah(No){Nh). Unlike the H2 molecules that desorb, the OH as well as H2O 
molecules (not included in the master equation above) are believed to stick to the grain and form 
an ice mantle. The parameters of the bare surface, used in this paper, are suitable only in the early 
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stages before the first monolayer of ice is formed. Beyond this stage one needs the activation energies 
for H diffusion on ice, which can be obtained from TPD experiments of hydrogen recombination 
on ice mantles. 

Note that in the equations above the coupling between different species is only through the 
expectation values (No) and (Ah). Therefore, one can simultaneously use the rate equations for 
some of the species and the master equation others according to the criteria discussed above. This 
is, in fact, an approximation in which the correlation between the probability distributions of 
different species such as Po(No) and Ph(Ah) is neglected. In cases where significant correlation is 
expected one can use a set of master equation for the joint probability distribution -Ph&o(^o> Ah). 
However, in this case the number of equations that are needed is N max (0) • N max (H), which may 
become impractical for systems with a larger number of species. In principle, such correlations 
could be significant in cases in which there are two species that react only with each other and thus 
the density of one specie would be strongly dependent on the availability of the other. However, in 
the interstellar medium such correlations are not expected to be strong since the reactive species 
(particularly atomic hydrogen) react with a large number of species. 

6. Summary 

We have introduced a new approach for the simulation of hydrogen recombination on micro- 
scopic dust grains in the interstellar medium. This approach is based on a set of master equation for 
the probabilities Ph(Nh), Ah = 0, 1, 2 . . . that there are Ah hydrogen atoms on the grain surface. 
Unlike the rate equations that provide a mean-field analysis, suitable for macroscopic surfaces, these 
rate equations provide an exact description of the recombination process on small grains taking 
into account the discrete nature of Ah as well as the fluctuations. The approach can be extended to 
more complex chemical reactions with multiple species and coupled to the densities of the reactants 
in the gas phase. 
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Figure Captions 



Fig. 1. — The recombination efficiency i], obtained from the master equation, for olivine grains 
of diameters d = 0.1/xm and d = O.Ol^m (x) under steady state conditions at a constant flux, 
as a function of the grain temperature. The flux is / = 0.18 x 10~ 8 MLs' 1 , which amounts 
to Fn = 10~ 3 and 10 -5 (atoms s _1 ) for the larger and smaller grain sizes, respectively. The 
recombination efficiency, obtained from the rate equations under similar conditions is also shown 
(solid line) showing a window of high efficiency between 7 - 9 K and a tail of decreasing efficiency 
above 9 K. We observe that in the center of the high efficiency peak there is good agreement 
between the master equation and rate equation results. However, for the smaller grain at higher 
temperatures, the master equation predicts significantly lower recombination efficiencies than the 
rate equations. These deviations persist at temperatures as high as 11 and 12 K but cannot be 
seen in the Figure due to the limited resolution. 

Fig. 2. — The expectation value for the number (Njj) of H atoms on an olivine grain as a function 
of the grain diameter under steady state conditions at T = 10K, as obtained from the master 
equation (O)- The solid line is simply a guide to the eye. 

Fig. 3. — The recombination efficiency r] for an olivine grain as a function of the grain size (•). at 
T = 10K and / = 0.18 x 10 -8 MLs^ 1 . The results predicted by the rate equations under similar 
conditions, which are independent of the grain size, are also shown (dashed line). It is observed 
that as the grain size decreases below d = 0.1/im the recombination efficiency quickly decreases. 

Fig. 4. — The distributions -Ph(-^h) on olivine grains of diameters d = 0.1 fim (O) an d d = 0.01 
/im (x), obtained from the master equation at T = 9K and / = 0.18 x 10~ 8 MLs^ 1 . 
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